#include<iostream>
//#include<pthread.h>
#include<immintrin.h>
#include<cmath>
#include<sys/time.h>

using namespace std;
#pragma comment(lib, "pthreadVC2.lib")

struct ThreadData {
	int start;
	int end;
	int k;
	int N;
	float** A;
};

//ÆÕÍšÏûÔªËã·š
void Lu(int N, float** A);
void print(int N, float** A);
void* parallelGaussianElimination(void* arg)
{
	ThreadData* data = (ThreadData*)arg;
	for (int i = data->start; i < data->end; i++) 
	{
		/*for (int j = data->k + 1; j < data->N; j++) 
		{
			data->A[i][j] = data->A[i][j] - data->A[i][data->k] * data->A[data->k][j];
		}*/
		////
		__m128 vaik;
		vaik = _mm_set1_ps(data->A[i][data->k]);
		int j;
		for (j = data->k + 1; j + 4 <= data->N; j += 4)
		{
			__m128 vakj, vaij, vx;
			vakj = _mm_loadu_ps(&data->A[data->k][j]);
			vaij = _mm_loadu_ps(&data->A[i][j]);
			vx = _mm_mul_ps(vakj, vaik);
			vaij = _mm_sub_ps(vaij, vx);
			_mm_storeu_ps(&data->A[i][j], vaij);
		}
		for (; j < data->N; j++)
		{
			data->A[i][j] = data->A[i][j] - data->A[i][data->k] * data->A[data->k][j];
		}
		////
		data->A[i][data->k] = 0;
	}
	pthread_exit(NULL);
}
void Lu_(int N, float** A) 
{
	pthread_t* threads = new pthread_t[N];
	for (int k = 0; k < N; k++) 
	{
		for (int j = k + 1; j < N; j++) 
		{
			A[k][j] = A[k][j] / A[k][k];
		}
		A[k][k] = 1.0;
		for (int i = k + 1; i < N; i++) 
		{
			ThreadData data;
			data.start = i;
			data.end = i + 1;
			data.k = k;
			data.N = N;
			data.A = A;
			pthread_create(&threads[i], NULL, parallelGaussianElimination, (void*)&data);
			pthread_join(threads[i], NULL);
			//std::cout << i << ":\n";
		}
	}
}

int main()
{
	struct timeval sign1, sign2, sign3;
	const int N = 10000;
	float** A = new float* [N];
	for (int i = 0; i < N; i += 1)
	{
		A[i] = new float[N];
		for (int j = 0; j < N; j += 1)
		{
			//cin >> A[i][j];
			A[i][j] = (i + 1) * (j + 1);
			if (i > j)
				A[i][j] = 0;
		}
	}
	//print(N, A);
	for (int i = 0; i < N; i += 1)
	{
		for (int j = 0; j < N; j += 1)
		{
			A[i][j] += A[i / 2][j];
		}
	}//³õÊŒ»¯ÂÒÅÅ
	//print(N, A);
	
	
	gettimeofday(&sign1, NULL);
	Lu_(N, A); cout << "Luing\n";
	gettimeofday(&sign2, NULL);
	
	Lu(N, A);
	gettimeofday(&sign3, NULL);
	cout << "p+s算法时间/us： " 
		<< (sign2.tv_sec*1000*1000+sign2.tv_usec)
			-(sign1.tv_sec*1000*1000+sign1.tv_usec);
	cout << " 平凡算法时间/us： " 
		<< (sign3.tv_sec*1000*1000+sign3.tv_usec)
			-(sign2.tv_sec*1000*1000+sign2.tv_usec);
	cout << "\n";
	
	//print(N, A);
	return 0;
}

void Lu(int N, float** A)//ÆÕÍšÏûÔªËã·š
{
	for (int k = 0; k < N; k++)
	{
		for (int j = k + 1; j < N; j++)
		{
			A[k][j] = A[k][j] / A[k][k];
		}//¹éÒ»»¯
		A[k][k] = 1.0;
		for (int i = k + 1; i < N; i++)
		{
			for (int j = k + 1; j < N; j++)
			{
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			}
			A[i][k] = 0;
			//cout << i << ":\n";
			//print(N, A);
		}
	}
}
void print(int N, float** A)
{
	for (int i = 0; i < N; i += 1)
	{
		for (int j = 0; j < N; j += 1)
		{
			cout << A[i][j] << ' ';
		}
		cout << '\n';
	}
}
